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Abstract 

Background: The Poisson-Boltzmann (PB) equation and its linear approximation have been widely used to 
describe biomolecular electrostatics. Generalized Born (GB) models offer a convenient computational approximation 
for the more fundamental approach based on the Poisson-Boltzmann equation, and allows estimation of pairwise 
contributions to electrostatic effects in the molecular context. 

Results: We have implemented in a single program most common analyses of the electrostatic properties of 
proteins. The program first computes generalized Born radii, via a surface integral and then it uses generalized 
Born radii (using a finite radius test particle) to perform electrostic analyses. In particular the ouput of the program 
entails, depending on user's requirement: 

1) the generalized Born radius of each atom; 

2) the electrostatic solvation free energy; 

3) the electrostatic forces on each atom (currently in a dvelopmental stage); 

4) the pH-dependent properties (total charge and pH-dependent free energy of folding in the pH range -2 to 18; 

5) the pKa of all ionizable groups; 

6) the electrostatic potential at the surface of the molecule; 

7) the electrostatic potential in a volume surrounding the molecule; 

Conclusions: Although at the expense of limited flexibility the program provides most common analyses with 
requirement of a single input file in PQR format. The results obtained are comparable to those obtained using 
state-of-the-art Poisson-Boltzmann solvers. A Linux executable with example input and output files is provided as 
supplementary material. 



Background 

Generalized Born models 

Electrostatic effects arising due to the interaction of 
solute charges among themselevs and with solvent and 
ion charges, are of utmost importance for biomolecular 
structure and function. Continuum methods based on 
the Poisson-Boltzmann equation have been widely used 
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for calculating electrostatic effects [1-5]. In the last dec- 
ades much interest has been devoted to developing fast 
approximations to the solution of the Poisson-Boltz- 
mann (PB) equation. 

Onufriev and coworkers have developed an analytical 
approximation to the exact potential inside and outside 
the low dielectric region of a sphere [6], that performs 
surprisingly well also for the complex shape of proteins 
[7-10] and is therefore more general than the general- 
ized Born (GB) models. 
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When only self- and interaction energies and forces 
are sought, the most widely used approach is based on 
generalized Born (GB) models. Recent reviews summar- 
ize the approach and highlight most interesting recent 
results [3,11-14]. 

Central to these models is the estimation of polariza- 
tion charge contributions to: i) the self-energy of each 
charge (embedded in the solute); ii) the interaction 
energy of each pair of charges. 

By equating the estimated reaction field \j™ act at the i th 
charge q t to the reaction field at a spherical ion in solu- 
tion with the same charge, each charge is assigned an 
effective radius which is called the Born radius (a,). 
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Once Born radii have been estimated the interaction 
energy between any two charges is computed as the 
sum of a direct Coulomb term (computed using the 
solute dielectric constant) and a solvation term which is, 

according to Still et al. [15], AG 50 / y = . AG so i Vi ij with: 
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This expression was found to be more accurate than 
the exact expression for two charges embedded in a 
conducting sphere, although the coefficient of 8.0 
instead of 4.0 in the exponential has been suggested [16] 
and other similar forms have been proposed [17,18]. 

Computation of Born radii 

Born radii could be in principle computed by solving a 
system of linear equations for the polarization charges 
at the boundaries of the solute volume [19-25]. Under 
the approximation that the ionic solution provides com- 
plete screening, amounting to the assumption that the 
surface behaves like a grounded conductor, polarization 
charges at the surface are such that the integral of their 
electric field at any outer surface point is exactly the 
opposite of the source charge field. This approximation 
amounts to setting s out - <*>. In a different context the 
same approximation has been proposed many years ago 
as the conducting surface model (COSMO) [26,27] and 
successfully used since then. Under this condition it is 
possible to obtain simple formulae for the generalized 
Born radius for a low dielectric region delimited by a 
sphere or a plane. The reaction field satisfies Laplace 
equation inside the surface and the solution may be 
obtained solving the equation using suitable basis 



functions with Dirichlet boundary conditions at discre- 
tized points on the surface [28]. 

These methods are however slower than approximate 
methods. The computation of Born radii is performed 
by volume or surface integrals. An approximation that 
has been widely used is the Coulomb field approxima- 
tion, where the formula for the electrostatic displace- 
ment for a uniform medium is applied in the inner and 
outer regions using the inner and outer dielectric con- 
stant, respectively. Born radius in this approximation 
does not depend on the inner and outer dielectric con- 
stant, but rather becomes a purely geometric quantity. 
The volume integral formulation of the Coulomb field 
approximation is turned into a surface integral formula- 
tion using the divergence theorem [29]. The Coulomb 
field approximation corresponds to neglecting the con- 
tribution to polarization at the surface due to polariza- 
tion charges and to considering therefore the electric 
field inside the surface as due only to the source charge: 
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Under the assumption of grounded conducting surface 
the outer field is zero and the polarization charge at the 
surface point r s is given by Gauss' theorem: 
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where n(?s) is the unit vector normal to the surface at 
point r s and pointing outwards. 

The reaction field at the source charge will be equal to 
the integral of the fields due to surface charges at the 
source charge, i. e.: 
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The reaction field may be used to compute the Born 
radius according to equation (1) with s out = °°. The 
Coulomb field approximation for the grounded conduct- 
ing surface has the advantage that the integral of the 
surface polarization charges will be always equal to the 
opposite of the charge inside the surface as for the exact 
solution. On the other hand it gives a very bad approxi- 
mation both of polarization charges, which depend only 
on the distance from the source charge and on the nor- 
mal to the surface, and reaction field energies even for 
simple systems. Although the Coulomb field approxima- 
tion is reasonable, Grycuk has pointed out that it is fun- 
damentally inaccurate for charges embedded in a sphere 
[30] by comparing both the self energy computed under 
the Coulomb field approximation and the pair interac- 
tion energy computed using the formula by Still and 
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coworkers [15], with the exact solution obtained by 
using the Kirkwood model [6]. Based on the analysis for 
charges embedded in a sphere an exact formula is pro- 
vided which may be recast as a function of a surface 
integral: 



jjreact _ _ I I 

4ne 0 s in \Js 



(r s - Tj) • n(r s ) dS \ 3 

wh-hw 6 4tt 



(6) 



from which the integral expression for the generalized 
Born radius follows: 
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The above expression which will be referred hereafter 
as GBR6 following Grycuk [30] and Tjong and Zhou 
[31,32], has been analysed in detail by Mongan et al [33]. 
Equation 7 was found to perform extremely well also for 
"very" non spherical shapes and in the context of biomo- 
lecular models. Occasional large differences with respect 
to Poisson-Boltzmann calculations were found for inner 
cavities and local concavities at the surface, i.e. in condi- 
tions where a continuum model is anyway questionable. 
That study concluded that with a correction for a small 
systematic error, the GBR6 model is a sufficiently accu- 
rate continuum electrostatic model [33]. 

Tjong and Zhou [31,32] used an analytical implemen- 
tation for the estimation of Born radii based on volume 
integrals (corresponding to the above formula (7) that 
uses a surface integral instead) and showed its superior 
accuracy compared to other existing methods for a set 
of 55 very different proteins. In their approach it is 
made clear that standard estimations of Born radii com- 
pute in fact only geometric properties, and they provide 
empirical formulae for the correct Born energy depend- 
ing on the inner and outer dielectric constants, ionic 
strength, total charge and number of atoms. 

Although successful and theoretically correct for a 
charge inside a grounded conducting sphere, there is no 
reason why equation (7) should provide good estimation 
of Born radii for complex shapes like those of proteins and 
biomolecules in general. In general Born radii could be 
estimated using integrals of the form (see the subsection 
Exact generalized Born radii for a conducting sphere): 
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In the absence of a theory which could be cast in a 
fast computational framework, fitting approaches have 
been successfully followed and tested on large sets of 
proteins showing excellent agreement with Poisson- 



Boltzmann calculation results. Romanov et al. [34] pro- 
posed that a linear combination of integrals I 3 to 7 6 , in 
the present notation, and a constant term could fit the 
self-polarization energy and thus be used to compute 
generalized Born radii (see also the discussion by 
Mongan et al. [33]). 

Applications of the generalized Born model 

Based on the correct estimation of Born radii, the solva- 
tion contribution to the interaction between any two 
charges may be computed by using equation (2). Deriva- 
tion of equation (2) with respect to atomic positions 
gives the electrostatic solvation forces acting on atoms. 
The implicit dependence of Born radii on atomic posi- 
tions makes computation of forces far from trivial 
[35-40] for the generalized Born model as well as for 
the parent Poisson-Boltzmann model where the bound- 
aries depend on atomic positions [41]. The possibility of 
computing energies and forces faster with respect to the 
reference Poisson-Boltzmann equation has made the 
Generalized Born model the choice of election for impli- 
cit solvent molecular dynamics simulations. Also, the 
computation of pairwise solvation energies allows for 
fast computation of pKa of multiple titrating groups as 
we discuss in the Methods section. 

Applications of generalized Born model (and other 
implicit solvent methods) have been reviewed elsewhere 
[3,11-14]. At variance with the reference Poisson-Boltz- 
mann model, the computation of electrostatic potential 
in space is outside the scope of the Generalized Born 
model where only interactions are considered thorugh 
equation 2. Here we use a finite radius test charge in 
order to define a potential within the frame of the gen- 
eralized Born model. 

Aim of this work 

In the present work we provide a software that, based 
on generalized Born radii, provides most common elec- 
trostatic analyses of proteins. The program first com- 
putes generalized Born radii, via surface integrals, based 
on different definitions and then it uses generalized 
Born radii (using a finite radius test particle, where 
needed) to perform electrostic analyses. In particular the 
ouput of the program entails, depending on user's 
requirement: 

1) the generalized Born radius of each atom; 

2) the solvation electrostatic free energy; 

3) the electrostatic forces on each atom (the theory of 
electrostatic forces based on surface integrals is devel- 
oped here and implementation is still at a developmen- 
tal stage); 

4) the pH-dependent properties (total charge and pH- 
dependent free energy of folding in the pH range -2 
to 18); 
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5) the pKa of all ionizable groups; 

6) the electrostatic potential at the surface of the 
molecule; 

7) the electrostatic potential in a volume surrounding 
the molecule. 

Results and Discussion 

Generalized Born radii 

The Born radii were computed from numerical surface 
integrals and the resulting self energies were compared 
with the reference ones obtained from the Poisson- 
Boltzmann equation for 1000 atoms randomly chosen in 
the test set of 55 proteins. Born radii and Poisson-Boltz- 
mann self-energies have been computed using the sol- 
vent accessible surface as dielectric boundary or the 
molecular surface computed using the program MSMS 
[42]. The surface points density used for MSMS compu- 
tations was 10 pts A" 2 . Table 1 reports the results con- 
cerning the self-energies obtained using different ways 
to estimate the generalized Born radii. It is seen that the 
GBR6 model performs very well and that the gain in 
accuracy when linear combinations of a constant term 
and self-energies corresponding to radii a n are used, 
with n ranging from 3 to 6 (LC5) and 3 to 10 (LC9), is 
limited. 

For best performance it is necessary to use a number 
of probe points per atom larger than 100. For less dense 
surface points occasional negative sign integrals are 
observed which need an adhoc treatment, e.g. resetting 
the generalized Born radius to a predetermined value. 
Although negative radii are not found for higher densi- 
ties of surface points, there are still rarely occurring 
large deviations between GB and PB self-energies, which 
could be due to differences in the computation of the 
boundary surface and its further treatment in the two 
approaches. As an illustration of the accuracy obtained, 
Figure 1 reports the plot of generalized Born self-ener- 
gies versus those computed using APBS. 



Table 1 Self energies 



model 


RMSD 


corr. coef. 


GBR6 SAS 


10.0 


0.995 


LC5 SAS 


3.5 


0.996 


LC9 SAS 


3.5 


0.996 


GBR6 MS 


57.6 


0.949 


LC5 MS 


23.9 


0.950 


LC9 MS 


22.7 


0.955 



Comparison of self-energies computed using generalized Born radii from 
surface integrals and using the Poisson-Boltzmann equation. For each 
computational model (SAS stands for solvent accessible surface, MS stands for 
molecular surface, see text) the average root mean square difference (RMSD) 
(kJ/(q mol)) and the correlation coefficient are reported. 
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Figure 1 GB vs. PB self energies. Generalized Born self-energies 
using solvent accessible surface integrals versus Poisson-Boltzmann 
self-energies computed using the program APBS. 



Electrostatic solvation free energy 

The electrostatic solvation energy computed using the 
solvent accessible surface integral (model GBR6) at each 
atom, computed as half of the product charge times 
reaction field, was compared with the corresponding 
solvation energy computed using the Poisson-Boltzmann 
equation. Compared to the solvation self-energy, the sol- 
vation energy depends also on other charges in the 
molecule and provides a test for the solvation effects 
computed according to equation 2. The average root 
mean square deviation from reference Poisson-Boltz- 
mann calculations is just 2.3 kJ/mol and the correlation 
coefficient is 0.9996. The accuracy of the GB model ver- 
sus the Poisson-Boltzmann model may be judged from 
Figure 2 which reports the solvation energy (computed 
as half the product of charge times reaction field) at 
each atom for the 93240 atoms of the 55 protein set. 

Electrostatic solvation forces 

Since solvation energies depend, according to equation 
2, on atomic positions explicitly and implicitly through 
Born radii, the computation of solvation electrostatic 
forces is quite complex, and strongly dependent on the 
interface model chosen [38-40] (see also for a general 
discussion [41]). 

In order to estimate how important are effects due to 
the dependence of Born radii on atomic positions, solva- 
tion electrostatic forces have been computed for each 
atom of the 55 proteins as the derivative of the solvation 
energy under the approximation that Born radii are con- 
stant. Albeit approximate this way of computing 
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Figure 2 GB vs. PB solvation energies. Generalized Born solvation 
energies at each atom using solvent accessible surface integrals 
versus the corresponding Poisson-Boltzmann energies computed 
using the program APBS. 



electrostatic forces preserve by definition zero total elec- 
trostatic force as expected for isotropic media and as 
found by the correct expression for the force [41]. Due 
to the different way of computing forces we did not 
attempt comparing ionic boundary, dielectric boundary 
and charge times electrostatic field components of the 
electrostatic force, but we rather compared the total sol- 
vation force. The results are not very accurate with the 
average root mean square deviation, with respect to the 
forces computed using APBS using the same para- 
meters, equal to 2.7 kJ/(Amol) compared to an average 
square root value of the force of 5.0 kJ/(Amol). The cor- 
relation coefficient is 0.84. 

The full computation of electrostatic solvation forces 
using surface integrals is described in the Methods sec- 
tion. As a demonstration of the accuracy of the calcula- 
tion 20 random point charges ranging between -1.0 and 
1.0 q have been embedded in a sphere of radius 5.0 A at 
a distance of at least 0.5 A from each other. Besides 
small differences due to treatment of the boundary, the 
agreement between force components computed using 
surface integrals and those computed solving the Pois- 
son-Boltzmann equation is very good, with a correlation 
coefficient of 0.9998 and a fitting slope of 0.983 (Figure 
3). Work is under way to implement force calculation in 
an efficient way for large systems like proteins. 

pH-dependent properties (total charge, pH-dependent 
free energy of folding, and pKa of ionizable groups 

The test set provided by Gunner and coworkers has been 
used to test the prediction of pKa's in proteins [43]. We 



< 50 




§ " 20 S()0 -150 -100 -50 0 50 
PB force components (kJ/mol A) 

Figure 3 GB vs. PB forces. Generalized Born force components at 
20 random points within a 5.0 A sphere with random charges, 
computed using solvent accessible surface integrals versus the 
corresponding Poisson-Boltzmann energies computed using the 
program APBS. 



compared the results obtained with the results obtained 
by running the program Propka2.0 [44-47]. The latter 
program is very fast and provides predictions whose 
accuracy is comparable to that of more computationally 
intensive programs. We chose this program as a refer- 
ence because it is widely used and it is seemingly the fast- 
est available with a very good accuracy. A wide range of 
programs and methods have been compared recently 
[48] and the reader is referred to that work for more 
extensive comparisons of existing softwares. 

The results are summarized in Table 2. It should be 
noted that the method used here was not extensively 
optimized to maximize correlation with experimental 
results or to minimize the root mean square deviation 
(RMSD) with respect to the experimental data. 

On the other hand some scaling of contributions is done 
and therefore large deviations from experimental results 
are not found. In this respect, Propka that uses ten adjus- 
table parameters and other parameters chosen for best 
performance [46] apparently does not prevent very large 
shifts to be predicted. As a consequence the global RMSD 
from experimental data is large and could thus easily be 
reduced. The execution time of Propka is in the range of 
seconds, while our program runs in minutes. No optimiza- 
tion or approximation has been implemented as yet. 

The results are somewhat better than those obtained 
by Propka (v. 2.0), although the figures reported here 
for the performance of Propka (v. 2.0) are worse than 
those reported in the original papers, because of a dif- 
ferent test set. We remark that the comparison reported 
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Table 2 Predicted vs. experimental pKa shifts 







BLUUES 


Propka v. 2.0 




RMSD 


corr. coef. 


mvl jU 


corr. coer. 


rOr 


n 79. 

u./ o 


0 42 


1.50 


0.27 


GLU 


0.90 


0.54 


1.75 


0.24 


HIS 


0.98 


0.74 


1.50 


0.59 


LYS 


0.79 


0.45 


0.79 


0.25 




1.78 


0.41 


2.12 


-0.03 


NTR 


0.54 


0.93 


0.96 


-1.00 


CTR 


1.01 


-0.14 


1.06 


0.36 


All 


1.00 


0.51 


1.59 


0.36 



Comparison of predicted vs. experimental pKa's. The RMSD is reported 
together with the correlation coefficient. We report the results obtained using 
Propka v.2.0. The latter program is very accurate, but occasional large shifts, 
not filtered nor treated ad-hoc here, result in large RMSD and smaller 
correlation coefficients. 



in Table 2 is dominated by outliers, that could be easily 
filtered or treated in an ad hoc manner, for Propka. 
Similar figures are found in independent tests [48-50]. 
An additional problem is that the data set includes mul- 
timeric (up to decameric) proteins and therefore the lat- 
ter contribute more than others to the figures reported 
in Table 2. It should be noted that this comparison is 
far from exhaustive as many other features could be 
chosen to judge a method's value. Extensive compari- 
sons have been performed by Stanton and Houk [48]. 
The purpose of the comparison performed here is to 
demonstrate that our program outputs results compar- 
able (or better on the dataset used here) than the results 
obtained by one of the best state-of-the-art method. In 
order to better judge the performance of the method, 
computed versus experimental shifts in pKa values, with 
respect to the reference model ones, are reported in Fig- 
ure 4. 

Surface electrostatic potential 

The evaluation of the potential in regions outside the 
molecule is beyond the scope of the generalized Born 
approach which focuses instead on self and interaction 
energies. This is at variance with approaches that approx- 
imate the potential computed using the Poisson-Boltz- 
mann equation. In particular Onufriev and coworkers 
[7,8] have found an approximation (not simply amount- 
ing to a truncation) based on Kirkwood's series expan- 
sion of the solution for a sphere [6]. Their analytical 
model performs surprisingly well for a large set of pro- 
teins and enables fast calculation of the potential at the 
surface and in the volume surrounding the molecule. 

The surface potential has been computed here, within 
the framework of generalized Born approach, as the 
energy of interaction of all atoms of the molecule with a 
unit test charge with a generalized Born radius equal to 
half the solvent probe radius, i.e. 0.7 A. This value was 




-4 -2 0 2 4 

Experimental pKa shift 

Figure 4 GB vs. PB pKa shifts. Computed versus experimental pKa 

shifts (compared to the reference model ones) for the dataset of 

Gunner and coworkers [43]. 
v ) 

found to reproduce well the potentials computed using 
APBS. Each surface point is assigned to the atom con- 
tacted by the solvent and for each solvent exposed atom 
the average surface potential is listed in the beta factor 
field in the output pdb file. The result is thus a readable 
list of atoms with the corresponding average surface 
potential if the atom is solvent exposed. This information 
may be used to display the potential using softwares like 
VMD [51]. As an example the potential computed using 
the Poisson-Boltzmann equation and using the method 
described above is shown in Figure 5 for the paired 
domain of the protein PAX6 and its cognate DNA. 
Exactly the same parameters are used in the two calcula- 
tion with the only difference that the surface in the Pois- 
son-Boltzmann calculation is generated using van der 
Waals radii inflated by 0.7 A. The boundary surfaces are 
however different in the two calculations and this value 
was found to compensate for the differences. 

For the same reason a point by point comparison of the 
potential at the surface is not completely appropriate, 
because boundaries are slightly different in the two calcu- 
lations. Instead of taking the average surface atomic 
potential, we compared the potential computed using the 
generalized Born model and a 0.7 A radius test charge at 
each surface point with the Poisson-Boltzmann potential 
computed by UHBD at the same point. Such comparison 
is reported in Figure 6 and, notwithstanding the differ- 
ences in methodologies, the correlation coefficient is 0.90 
and the average error is just +0.08 kcal/(mol q). The root 
mean square deviation is 0.43 kcal/(mol q) with the lar- 
gest contributions to this figure due to outliers by several 
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Figure 5 GB vs. PB surface potential (average over atomic surface). Computed surface potential espressed in kcal/(mol q) for the paired 
domain of the protein PAX6 and its cognate DNA (PDB structure id.: 6PAX). The reference Poisson-Boltzmann potential (lower row) and the 
average generalized Born potential at atomic surface (upper row) are displayed with colors ranging from red (-3.0 kcal/(mol q)) to blue (3.0 kcal/ 
(mol q)) for DNA and from red (-2.4 kcal/(mol q)) to blue (2.4 kcal/(mol q)) for the protein. 



standard deviations. Outliers are found in inner cavities, 
but also at concavities at the exposed surface. In some 
cases the errors are due to the differences in boundary 
definition. The error distribution is however similar to 
that reported by Onufriev and coworkers [8]. 

Electrostatic potential in a volume surrounding the 
molecule 

Similar to the surface potential the electrostatic poten- 
tial for a unit test charge with generalized Born radius 



equal to half of the solvent probe radius may be com- 
puted at nodes of a grid enclosing the molecule and 
output in the DX format (the DX file format is 
described at [52]) readable by programs like VMD [51]. 
The grid is chosen based on as the minimum box entail- 
ing the molecule plus three Debye length on each side. 
If the Debye length is more than 10 A, 30 A are added 
at the minimum box on each side. The space is then 
divided in 97 x 97 x 97 grid points and the potential is 
computed at each point. The method is thus slower 
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PB surface potential 

Figure 6 GB vs. PB surface potential. Surface potential, espressed 

in kcal/(mol q), computed using the generalized Born model, versus 

the reference Poisson-Boltzmann potential at the same surface point, 
v ) 



compared to solving the finite difference equation for 
the potential. At variance with standard options with 
other software, here the potential inside the molecule is 
set to 0, which removes isopotential curves inside the 
molecules (uninformative because arising from strong 
varying local fields) and helps with the visualization of 
the structure of the molecule together with outer isopo- 
tential curves. Example input files are given in the sup- 
plementary material. The comparison with the 
analogous analysis using the program APBS shows less 
smooth isopotential surfaces (Figure 7), a fact which 
might be due to the finite radius of the unit charge test 
particle. Note that the output potential is in kj/(q mol) 
compared to widely used kT/(q mol) or kcal/(q mol). 

Conclusions 

A program for the analysis of the electrostatic properties 
of proteins based on the surface integral computation of 
generalized Born radii has been presented. Further work 
will be devoted to improve the efficiency of calculations 
in particular for what concerns electrostatic solvation 
forces. The program, together with examples is given as 
supplementary material (Additional file 1). 



temperature 298.15 K, ionic strength 0.150 M. In APBS 
the mesh of the grid was 0.25 A which resulted in very 
large grid size (257 3 points). In UHBD, used for surface 
potential calculations, a focusing procedure was used 
with the final mesh size of 0.8 A. 

The boundary surface was defined in the two pro- 
grams in order to match as close as possible the defini- 
tion used in our program, i.e. the molecular surface for 
atoms with radii inflated by the radius of the solvent. 
The surface point density was 10 A" 2 . 

Except where noted the same temperature, dielectric 
constant and ionic strength were used in our program. 

Solvent accessible surface 

Solvent accessible surface points and surface normal 
vectors have been generated by considering the van der 
Waals sphere of each atom inflated by the radius of the 
solvent (e.g., for an atom of radius 1.9 A and a solvent 
radius of 1.4 A the inflated radius was 3. 3 A). 

Points on the inflated van der Waals sphere and sur- 
face normal vectors are precomputed for a set of 
inflated radii (r) and centered at atom positions. The 
coordinates of the i-th point (i ranging from 1 to n pt ) 
on the sphere are generated according to the following 
rule based on the golden section (due to Anton Sher- 
wood, see http://www.cgafaq.info/wiki/Evenly_distribu- 
ted_points_on_sphere and links thereof): 

Zi = r x 1 -2 2 - 

\ n pt ) 

Xi = tJt 2 — zf x cos((3 — V5)yr(i — 1)) 
Yi = Jr 2 - zf x sin((3 - V5)jr{i - 1)) 

n pt was 200 corresponding for an atom with radius 
1.9 A to a density of 4.4 pts A" 2 . 

Atoms are mapped on a grid and the position of each 
surface point is compared only with the position of the 
list of atoms associated with neighboring grid points. 
This procedure was found to be efficient and robust 
without any failure. 

Each surface point within the inflated van der Waals 
volume of a different atom is considered buried. All sur- 
face integrals have been performed as the corresponding 
finite sums over exposed (non buried) surface points. 



Methods 

Reference APBS and UHBD calculations and parameters 

Reference electrostatic calculations were performed with 
the programs APBS [53] and UHBD [54,55]. Except 
where noted, standard parameters were used: inner 
dielectric constant 1.0, outer dielectric constant 78.54, 



Generalized Born radii 

Reference Born radii have been obtained from Poisson- 
Boltzmann self-polarization energies computed using 
APBS [53] according to eq. 1. 

Different Born radii have been obtained from surface 
integrals as: 
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Figure 7 GB vs PB isopotential curves. Isopotential curves at 0.3 kcal/(mol.q). Upper row: potential obtained using the GBR6 model; lower row: 
potential obtained by solving the Poisson-Boltzmann equation. 




where I n is defined by equation 8, or by fitting the 
reference solvation self energy by linear combinations of 
the solvation energies corresponding to different oc n 's (n 
ranging here from 3 to 6) as done by Romanov et al. 
[34]. Fitting was necessary because the published coeffi- 
cients did not provide good results, as a result of the 
different procedures used to compute surfaces and sol- 
vation free energies. 

An equation for the generalized Born radius 

The reaction field obeys Laplace equation in the solute 
region. When the reaction field is expressed via a gener- 
alized Born equation like (2), we may conveniently study 
what are the requirements imposed by Laplace equation 
on Born radii. In particular we consider the reaction 
field due to a source charge at r { and consider the reac- 
tion field at a point rp Lly, in the limit r ; Under 
this limit the reaction field may be approximated as: 



where the constant K depends on the coefficient used 
in the exponential and is 0.75 for the commonly 
adopted Still formula. 



Imposing that V 2 i/,y = 0, where all derivatives are with 
respect to coordinates of rp and letting r ; tend to fp after 
some straightforward manipulation, results in the fol- 
lowing equation for the Born radius a t \ 

aV 2 a - ~(Va) • (Va) + 6K = 0 (10) 

For the conducting grounded surface approximation 
the boundary condition is that the Born radius 
approaches 0 as the boundary surface is approached. It 
is easy to check that this equation is satisfied for the 
sphere and the plane where K is equal 1. 

For two parallel conducting infinite planes (approxi- 
mated here by two very large circular plates) the solu- 
tion may be compared with that obtained by expanding 
the exact solution in Bessel functions. We set the inter- 
plates distance along the z-axis equal to 1.0 and the 
radius R of the plates equal 1000.0. Consider x as the 
distance of the source charge from one of the plates on 
the axis of symmetry. The reaction field between the 
plates may be expressed as a sum of functions which 
satisfy the Laplace's equation in cylindrical coordinates: 

UreactiPrZ) = J> fe + ^ + C~ <T te )/o(fep) (U) 

k 

where J 0 (x) is the Bessel function of order zero. We 
may consider only order 0 functions due to cylindrical 
symmetry. The coefficients are obtained by imposing 
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/ 



that at the boundary the potential (U react + U source ) is 
equal to zero and using the orthogonality relation: 

xJo(kx)Jo(k f x) = \&{k — k'\ The parameters k are 
k 

chosen such that J 0 {kR) = 0. 

The solution of equation (10) is obtained by bisection, 
guessing first a value at the midpoint of the axis of sym- 
metry and integrating the equation using an adaptive 
Runge-Kutta fourth order method [56] and requiring 
that the value of the Born radius at the boundary be 0. 

The two solutions reasonably agree, within numerical 
accuracy, and are compared in Figure 8. Solving Equa- 
tion (10) for proteins could represent an alternative to 
other well established methods to compute generalized 
Born radii. Obviously all the above derivation rests on 
the assumption that Still's formula (or alike) provides a 
good approximation of the exact solution to the Poisson 
(or Poisson-Boltzmann) equation for non-homogeneous 
media. 

Electrostatic solvation energies in ionic solutions 

A merit of generalized Born models is to separate Cou- 
lombic interactions from reaction field effects. The gen- 
eralized Born radii have been used to compute 
interactions in ionic solutions assuming that charges 
whose distance is much larger than their Born radii 
interact through a Debye-Huckel potential and that the 
self-polarization energy can be computed according to 
the Debye-Huckel law [57]. The expression used in the 
present work for the reaction potential due to a unit 



source charge i at the site ; for an ionic solution of 
dielectric constant s out and Debye constant k D is: 



JjSOlV 



1 


1 














47T£ 0 \ 





oi\e-. — - 



(12) 



which is slightly different from that reported by Srini- 
vasan et al. [57] in order to include dielectric constant 
different from 1. When we fit the self-polarization 
energy obtained by the above equation versus that 
obtained by equation (22) in the work Tjong and Zhou 
[31], using a zero intercept line and assuming a physio- 
logical ionic strength of 0.150 M we find a correlation 
coefficient of 0.999 and a slope of 0.999. The total solva- 
tion energy is obtained by summing contributions from 
all pairs of atoms and solvation self-energies: 



(13) 



= 1,N 



The total energy of the molecule, which does not 
include self-energies in the inner dielectric medium is: 

Gun = ^ E MAi = ^ E m ( 4n£ J r . + K) (14) 
z i,i=i,N / _i=i m \4ne 0 e m r t j / 



ij'=l,N 



where the functions (p it j are the Green's functions at 
points i and j but do not include the self potential of 
each charge in the inner dielectric medium and S if j is 
Kroneckers delta. 




x 



Figure 8 GB radius from Equation (10) vs. numerical solution. 

Computed generalized Born radius for a charge between two 
infinite planes depending on the distance x from one of the planes. 
The generalized Born radius and the distance are normalized by the 
plane to plane distance. Solid line refers to the solution obtained 
from equation (10) and dashed line refers to the expansion (11) of 
the solution for a finite system. 



Electrostatic solvation forces 

The electrostatic solvation forces are obtained by differ- 
entiation of G soiv with respect to atomic coordinates. 
The implicit dependence of Born radii on atomic coor- 
dinates through a surface integral makes the derivation 
rather tedious. One approximation is to consider that 
effects from the variation of Born radii with atomic 
coordinates are smaller than the explicit variation of 
energies due to the change in interatomic distances. 
This is in general not the case. Just consider the force 
arising from variation of the self-energy with atomic 
coordinates. For a point charge q in a conducting sphere 
of radius R at a distance a from its center the solvation 
force is pointing away from the center of the sphere and 
its magnitude is: 

q 2 aR 

F = — 1 - j 15 

For a unit charge placed at 7.5 A from the center of a 
globular protein of radius 10 A, assuming an inner 
dielectric constant 4.0, the resulting force is 14 kJ/(mol 
A). An opposite force of the same magnitude is applied 
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on the center of the sphere defining the boundaries. 
This simple example shows that forces arising from var- 
iations of Born radii are not negligible. 

The derivative of electrostatic potential implies deriva- 
tives of Born radii with respect to atomic coordinates 
that may be computed as derivatives of a function of a 
surface integral: 



dxk 



dxk 



(16) 



The derivative of I l 5 is not straightforward as the 
domain of integration depends on atomic coordinates. 

The subtleties of derivation of integrals have been 
described by Flanders [58]. We use here a form valid for 
the derivation of a flux through a surface that does not 
entail the derivative of the boundary. The total deriva- 
tive with respect to a variable t (typically time, here 
atomic coordinates) of the flux of a vector field p 
through a surface S is expressed as a surface integral 
[59]: 



L 



dF 



+ (V • F)v(r s ) + V x (F x v(?s)) 



■n(r s )dS (17) 



where 



v(rs) = 



dfs 
dt 



The surface point (for the solvent accessible surface) is 
dependent on the contacting atom (say the j th atom): 

r s = fj + n{js){r vd w,j + r p ) (18) 

where r vdW j is the van der Waals radius of atom ; and 
r p is the radius of the solvent. The derivative of r s with 
respect to the k th coordinate of atom / takes a very sim- 
ple form: 



dfs 
dxi r k 



(19) 



where eu is the versor of the k th axis. Further deriva- 
tions of r s will be zero. Taking into account the latter 
statement and expanding the third term in the right- 
hand side of equation 17 a simpler form of the same 
equation is found: 



(rs) • V)F 



n(r s )dSi20) 



The left-hand side of equation 20 is directly related to 
equation 16 with: 



1 



4tt I In 



We define a function is (js) that assigns to each sur- 
face element the index of the atom defining the bound- 
ary and consider the derivative of I l 5 with respect to the 
lt h coordinate of atom /: 



dx ti 



i-/.K(iifcsF)-*™(fi'«)(il^i.)]-™ s < 2 W 



■M 



nk{r s ) 6 {fi-f s )-n{f s ) 



47r||r I --r s || 6 4tt |ft-r s || 8 



+ implicit surface-dependent terms 



i-h^dS (22) 



(23) 



Note that in the above equations the operator y oper- 
ates only on r s coordinates, i.e. at the point where the 
integrand function is evaluated. 

Let us consider now a simplified form (exact for a 
conducting sphere) of the Green's function for the elec- 
trostatic solvation energy 



1 



(24) 



and calculate explicitly the solvation force along the 
h th coordinate of atom ;: 



fik 



dG solv d(\Y*. } W}Ui4) 



1 



dxik 



dxik 



(Xj,k - Xj,k) 

qrl+aiaj) 



(25) 



(26) 



12 4^(^17^) 



o dlr o dl r 

5 '% +ff i&) (27) 



where the apices i and ; in I 5 means that the integral 
involves the vector from the surface points to the coor- 
dinates of atoms i and ; respectively. The computation 
of forces via surface integral according to the above 
equations is not practical, unless a cutoff on the distance 
between pairs of atoms and surface points is used. 
N(N + I) 

Indeed for each of the — pairs of atoms ; and i, 

where Af is the number of atoms, surface integrals must 
be computed (see equation (21)). This is the result of 
the brute force application of our surface integral defini- 
tion of Born radius. Work is under way in our labora- 
tory to make this computation efficient. 
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pH-dependent properties (total charge and pH-dependent 
free energy of folding in the pH range -2 to 18) 

The method described by Antosiewicz et al. [60] is 
implemented here with some modifications. PDB entries 
corresponding to proteins for which data are available in 
the pKa database made available by Gunner and cowor- 
kers [43] were prepared with pdb2pqr [61]. 

The method implemented here is a single conformer 
method, at variance with other methods like the QM/ 
MM/TI method by Simonson et al. [62], MCCE [43] 
and the constant pH molecular dynamics method by 
McCammon and coworkers [63]. These methods are 
expected to be more accurate, at the expense of large 
computational time. 

Approaches similar to that implemented here have 
been used by Onufriev and coworkers [64,65] where a 
significant speedup for larger proteins, with respect to 
our implementation, is due to clustering of interacting 
sites. 

Significant improvements over continuum electro- 
statics (and more in general molecular mechanics) 
methods have been achieved by parametrization of 
interactions [49,50,66-69]. 

Other approaches which use a continuum method 
have been proposed recently which achieve very good 
agreement with experimental values, by using a pentape- 
tide reference model for the titratable group, and by 
optimizing the hydrogen positions [70]. 

Improvements similar to those mentioned in the 
above paragraphs will be implemented in the future in 
our program. 

We summarize here the theory that is discussed at 
length by Antosiewicz et al. [60]. 
Energy of ionization 

The free energy of ionization of an ionizable group in an 
unfolded polypeptide at a given pH is: 

AG = \og{10)RTz{pK a - pH) (28) 

where z is the charge of the group upon ionization, R 
is the gas constant, T is the temperature (298.15 K in 
the present calculations). 

It is assumed that pKa of protein ionizable groups in 
an unfolded protein are the same as those of model 
compound. We assume here that the the pKa are the 
following: 4.5 for Glu, 3.8 for Asp, 6.5 for His, 12.5 for 
Arg, 10.5 for Lys, 9.0 for Tyr, 8.0 for Cys, 8.0 for the N- 
terminal ammine and 3.2 for the C-terminal carboxyl. 

The protein environment may change the free energy 
of ionization. If we assume that the main effects are 
electrostatic in nature the free energy of ionization for 
groups in a folded proteins is: 

AG = log(10)m f (p^ - P H) (29) 



+ \ Ota + Zf )ta + z i) - MjWij (30) 

ij 

~\ Ota + Zl ')ta + z i) - MjWij (31) 

ij 

where, for simplicity of notation, the sum runs on all 
atoms of the protein. q t are the atomic charges in the 
neutral amino acids and z t is 0 for most atoms and 1 or 
-1 for those representing basic and acid ionized groups, 
repectively. Superscript M stands for the model com- 
pound representing the unfolded state and N stands for 
the natively folded protein. The model compound of 
ionizable residues is obtained by excising that residue 
from the protein. The Green's functions q> it j are given by 
equation 14. The equation can be rearranged as follows: 

AG = £iog(io)ra^(p^ (32) 

i i ij i>j 

where each term describes a contribution to the free 
energy of ionization: the first term is the free energy of 
ionization in the model compound, the second term is 
the difference in ionization self-energy, the third term 
describes the difference in interaction of ionization 
charges with partial charges in the unionized protein 
and model compound and the fourth term describes the 
interaction among ionization charges. Following Anto- 
siewicz et al. the inner dielectric constant is set to 20.0. 
A wide range of inner dielectric constants has been used 
for pKa calculations. When a single dielectric constant 
has been used it has been chosen mostly larger than 1, 
e.g. 4 in MCCE [43], 8 in EGAD [68] and 11 in GB/ 
IMC [70]. A large value of the dielectric constant 
accounts in an empirical way for the missing degrees of 
freedom implied by a single conformer method (see the 
discussion by Schutz and Warshel [71]). 
Monte Carlo simulation of the ionization state 
The protonation state of the ionizable sidechains is 
changed according to a Monte Carlo procedure. Proto- 
nation or deprotonation are simulated by adding or sub- 
tracting a unit charge to an atom representing the 
ionizable group. The following atoms have been consid- 
ered as representative for each ionizable group: CG for 
Asp, CD for Glu, NZ for Lys, CE for Arg, SG for Cys, 
OH for Tyr, N for the N-terminal ammine and C for 
the C-terminal carboxyl. For histidine protonation may 
occur alternatively at the atom NE2 and ND1. 

Monte Carlo simulations are performed at 0.5 pH 
intervals between -2.0 and 18.0 and at each pH value 
average ionization values and average components of 
ionization free energy are computed. At each pH value 
the number of equilibration steps is 100 times the num- 
ber of ionizable sites and the number of Monte Carlo 
steps is 1000 times the number of ionizable sites. 
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Based on Monte Carlo simulation the output of the 
program includes: i) the titration curve for each ioniz- 
able site; ii) the list of pKa values for each ionizable site; 
iii) the charge state of the protein; iv) the pH-dependent 
component of the free energy of folding. Compared to 
the scheme of Antosiewicz et al. we are using the sol- 
vent accessible surface (as detailed in the subsection 
"Solvent accessible surface") instead of the molecular 
surface because results are more stable. 

Heuristic corrections (detailed hereafter) to the 
scheme of Antosewicz et al [60] have been done mainly 
to remove large contributions which are typically over- 
estimated due to neglection of molecular flexibility. For 
instance for surface residues large unfavorable interac- 
tions may be accomodated by conformational relaxation. 
Also large charge-charge interactions are likely to lead 
to partial unfolding that relieves the strong unfavorable 
energy. These considerations are consistent with the 
idea of using different dielectric constants for buried 
and solvent exposed residues, as suggested many years 
ago by Demchuck and Wade [72]. 
Scaling self energies 

Here we consider as solvent exposed all those sites that 
have generalized Born radius smaller than 4.0 A and 
buried those that have generalized Born larger than 7.0 
A. Self-energy interactions (second right-hand term in 
Equation 32) are multiplied by a factor 2.0 for buried 
sites and by 0.25 for exposed residues. For all intermedi- 
ate situations the multiplying factor is a linear combina- 
tion of the two extreme values. 
Scaling background interactions 

For exposed residues unfavorable background interac- 
tions (third right-hand term in Equation 32) are 
weighted by the empirical function exp(— ^). The fac- 
tor 2 in the Boltzmann-like function is purely empirical 
to downweigh a contribution that is likely to be relaxed 
by protein flexibility. 

For each ionizable site the pKa is obtained as the mid- 
point of the titration curve. The value is further cor- 
rected by separating the shift in pKa due to the 
desolvation self-energy {dpK s ^ )> background interactions 
(dpKa g ) an d site-charge - site-charge interactions (dpK 1 *). 
Final scaling of contributions 

The site-charge - site-charge interactions shift is scaled 

dpK il x 3 0 

by the function dpK" = a . — which sets 3 pKa 

units as the maximum contribution arising from interac- 
tions between titratable sites, consistent with the idea 
that large unfavorable interactions lead to partial or glo- 
bal unfolding. Similarly, background interaction shifts 

are scaled by the function dpKa g = which 

\dpK b a g \ + 5.0 



sets 5 pKa units as the maximum contribution due to 
background interactions, in order to avoid few very 
large shifts. 

Test sets and conditions 

The 55 proteins selected by Tjong and Zhou [31] for test- 
ing the GBR6 model have been used here. Atoms of the 
PDB structures have been assigned charges and radii taken 
from the CHARMM forcefield, using the program 
PDB2PQR [61,73]. The PQR file format is described at 
[52]. The radius of hydrogen atoms was reset to 1.0 A in 
order to avoid numerical inaccuracies in the analysis linked 
to the small radius of polar hydrogens in the CHARMM 
forcefield. The surface of the molecule was defined alterna- 
tively as the surface accessible to the center of a solvent 
probe sphere with radius 1.5 A or as the surface accessible 
by the surface of a solvent probe sphere with radius 1.5 A, 
except where noted. The first type of surface, here referred 
to as solvent accessible surface, is computed by default by 
the program while the second type of surface, here referred 
to as molecular surface, was computed using the program 
MSMS [42] and given as input to the program in the form 
of a list of vertices and normal vectors. 

For the computation of self energies 1000 sites were 
randomly chosen in the 55 proteins providing an 
unbiased sample of different environments. 

The inner dielectric constant was assumed to be 1.0, 
the outer dielectric constant was assumed to be 78.54, 
the temperature is 298.15 K, ionic strength 0.150 M, 
except where noted. For the computation of pKa of 
ionizable sites in proteins we used the database devel- 
oped by Gunner and coworkers [43] available at [74]. 
The test set includes pKa for structures with PDB id: 
la2p, la6k, lbeg, lbf4, lbhc, lbus, lbvi, lbvv, lcdc, 
lcoa, lcvo, lde3, ldg9, ldwr, legf, lgbl, lgoa, lh4g, 
lig5, ligc, lkf3, llni and Use. 

Exact generalized Born radii for a conducting sphere 

We analyze the simple planar and spherical systems 
which are treatable analytically with closed formulae in 
the hypothesis of grounded conducting surface. We con- 
sider a hollow conducting grounded sphere of radius R 
filled with a medium with dielectric constant s in . The 
reaction field for a point charge in the sphere may be 
easily obtained by the image charge method. The reac- 
tion field due to a source charge displaced by a with 
respect to the center of the sphere is computed as the 
field due to an image charge q f = —cj— at a distance 
from the center — along the direction from the center 
to the source charge: 

jjreact, sphere _ tfi ^ (33) 

1 AnsoSi n R 2 - a 2 



Fogolari et al. BMC Bioinformatics 2012, 13(Suppl 4):S18 
http://www.biomedcentral.com/1471-2105/13/S4/S18 



Page 14 of 16 



which corresponds to a Born radius: 

R 2 -a 2 



R 



(34) 



In the limit of R — > °° and a — > R the above expres- 
sion reduces to: 



react,plane 



1 



AnsQSi n 2d 



(35) 



where d is distance of the charge from the delimiting 
plane. The polarization charge contribution to the 
potential due to a source charge j at a point i inside the 
sphere is given by: 



Uij = - 



47tS 0 Si„Jjf j +aiOtj 



(36) 



where a t is the generalized Born radius at point i. 

As noted by Grycuk [30] for a charge placed at a dis- 
tance a from the center of a sphere with radius R, the 
Born radius computed according to equation (5) results in: 



1 



2^R 2 



1 R + a 

+ — lo B7; — :) 



la 



R-a J 



(37) 



This provides the exact Born radius in the limit a — > 
0, i.e. when the charge is placed at the center of the 
sphere, but underestimates Born radius in the limit R — > 
oo and a — > R (i. e. the plane) by a factor 2. It would be 
desirable to find a surface integral that recovers the cor- 
rect expression at least for charges embedded in a 
sphere. The same surface integral should have a form 
similar to that of equation (5) in order to ensure that no 
problems arise from complex shape surfaces, e. g. like 
those of biomolecules. A convenient correction may be 
sought considering other integrals like in equation (5) 
with increasing powers of the distance ||r f —r&\ \ In par- 
ticular we may define the following integrals: 



h 



I 4* 



(r s - n) • n(r s ) 



s 4tt \\rs-n\ 



I 



1 (r s - r,-) ■ n(r s ) 
' s 4jt Mrs — nil 5 



1 (js - n) ■ n{f s ) 
4ji \\r s -ri\\ n+1 



dS 



dS 



dS 



(38) 



(39) 



(40) 



For the sphere these integrals have closed forms that 
can be obtained using the formal integrator of Mathe- 
matica [75]. If the position of the source charge is dis- 
placed from the center of the sphere by a quantity a 
and R is the radius of the sphere the integrals have the 
following values: 



1 R 1 R+a 
h = 2 { W^ + 2a X ° g R-^ 



h = 



h = 



h = 



h = 



4a 2 



6{R 2 - a 2 ) 2 R 2 ~ a 2 
R 3 

(R 2 - a 2 ) 3 

24a 4 + 40a 2 {R 2 - a 2 ) + 15{R 2 - a 2 ) 2 
15{R 2 -a 2 ) 4 

16Ra 4 + 22a 2 {R 2 - a 2 ) + 6R(R 2 - a 2 ) 2 
6(R 2 - a 2 ) 5 



(41) 



(42) 



(43) 



(44) 



(45) 



The integral I 5 is interesting because, as pointed out 
by Grycuk [30] it can be easily inverted to get the cor- 
rect Born radius for both the plane and the sphere: 



1 

1 \ 3 



(46) 



The work of Mongan et al. [33] reports a compact 
form for the Born radii corresponding to the integrals 
above expressed in term of GBR6 multiplied by a cor- 
rection factor. 

Additional material 



Additional file 1: Executable and examples. Zipped folder with the 
executable compiled for Linux x86_64 and 80386 machines and 
including example input and output files. 
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PB: Poisson-Boltzmann; GB: Generalized Born; COSMO: Conducting Surface 
Model; RMSD: Root mean square deviation. 
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